Signatures of Spin and Charge Energy Scales 
in the Local Moment and Specific Heat 
of the Two— Dimensional Hubbard Model 
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Local moment formation driven by tlie on-site repulsion U is one of the most fundamental features 
in the Hubbard model. At the simplest level, the temperature dependence of the local moment is 
expected to have a single structure at T ^ U, reflecting the suppression of the double occupancy. 
In this paper we show new low temperature Quantum Monte Carlo data which emphasize that the 
local moment also has a signature at a lower energy scale which previously had been thought to 
characterize only the temperatures below which moments on different sites begin to correlate locally. 
We discuss implications of these results for the structure of the specific heat, and connections to 
quasiparticle resonance and pseudogap formation in the density of states. 
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I. INTRODUCTION 

The finite temperature properties of the two dimen- 
sional Hubbard model have beep extensively studied both 
analytically and numerically.El Quantum Monte Carlo 
(QMC) is especially effective at half-filling, where there 
is no sign problem. From calculations of the magnetic 
structure factor, susceptibility, compressibility, density 
of states, and the electron self-energy, a clear picture 
has emerged concerning the nature of the short and long 
range magnetic order, the Mott gap, and the quasiparti- 
cle dispersion. 

While the specific heat C{T) has been computed by a 
number of groups i^i one-dimension, principally by Bethe 
Ansatz technique^, there have been few QMC stud- 
ies of CXJ^ for the two and three dimensional Hubbard 
models.H~El The behavior for large U is well understood 
and one expects, as in the one dimensional case, a two 
peak structure in C'(T), with a broad high temperature 
peak at T ^ U associated with "charge fluctuations" , 
and a narrower peak at lower temperatures associated 
with "spin fluctuations" . 

This two peak structure of the Hubbard model can 
be understood from a strong coupling viewpoint as fol- 
lows: At temperatures which exceed the on-site repul- 
sion, T > U, the up and down electrons are decoupled, 
{n-ini) = {ni){ni) = j, at half-fiUing, and the local mo- 
ment (ml) = ((n| - n|)2) = {n^ + - 2n^ni) = 1 - 
{n^ni) = 0.5, its uncorrelated value. At the temperature 
scale T ^ U, double occupancy begins to be suppressed, 
(n|n|) — > 0, and 1. Since the potential energy 



in the Hubbard model is just P = Un^ni — ^(1 — 
this growth in the local moment is synonymous with a de- 
crease in the potential energy as T decreases, and hence 
a peak in C(T) at T ~ [/. Once moments are formed, 
the half-filled Hubbard model maps onto the spin-1/2 
antiferromagnetic Heisenberg model, with exchange con- 
stant J = At^ /U , whose specific heat has a peak at T J 
associated with magnetic ordering. In sum, one expects 
C{T) for the strong coupling Hubbard model at half- 
filling to have a "charge peak" at T ^ U and a "spin 
peak" at T J. This strong coupling argument further 
suggests that the temperature derivatives of the poten- 
tial and kinetic energies are associated with the high and 
low temperature specific heat peaks respectively. 

On the other hand, the behavior of the specific heat for 
small U in the two dimensional Hubbard model is still 
unclear. In particular, it is not known whether a two 
peak structure is present for all values of the interaction 
or if these coalesce into a single peak at small U . That 
the two peaks might merge is suggested by the fact that 
the charge and spin energy scales, U and /U, approach 
each other as U is decreased. However, the situation is 
not so straightforward, because the strong coupling form 
for the energy scale for SWi ordering crosses over to a 
weak coupling expressionO'D t exp(— 27r-\/t/?7), which at 
small U is still well-separated from the energy scales t 
and U. 

The QMC results reported in this work are consis- 
tent with such weak coupling behavior at small U . That 
is, one of our key findings is that C{T) shows two dis- 
tinct peaks which persist to couplings an order of mag- 
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nitude less than the noninteracting bandwidth. It may 
be that the half-fiUed two-dimensional Hubbard model is 
unique in this respect, since, as we review below, stud- 
ies in one dimension and within DMFT show a merging 
of the peaks at values of U roughly equal to the band- 
width. Indeed, the two-dimensional half-filled Hubbard 
model on a square lattice has unusual nesting features in 
its noninteracting band structure, as well as a logarith- 
mic divergence of its density of states at the Fermi level 
which have previously been noted to enhance antiferro- 
magnetism anomalously.Q 

Whether this is the case or not, there is as yet no com- 
pelling evidence of the appearence of a weak coupling 
energy scale texp{—2TTt/U) in the specific heat in other 
dimensions. Jn one dimension, exact pdiagonalization in 
small chainscl and QMC calculations^ suggest the two 
peaks merge, but disagree with respect to the interaction 
strength at which this occurs. Exact diagonalization is 
limited to chains of very modest extent, and finite size ef- 
fects tend to be large for small U. Meanwhile, QMC work 
did not reach low enough temperatures to resolve the two 
peaks even for large values of U where they almost cer- 
tainly both exist. Bethe-Ansatz calculations Jaelp clarify 
this issue, but focus on the large U limitO^EJ. Despite 
these various caveats, the consensus of these approaches 
is that the spin and charge peaks merge at U/t « 4, the 
one-dimenskmal bandwidth. Quantum transfer matrix 
calculationgEj also show the merging of the two peaks at 
U/t ^ 4. 

Coalescence of the specific heat peaks has also been 
seen in "Dynamical Mean Field Theory" (DMFT) 
studies the system in the limit of high dimension.EfrE£l 
There, the Hubbard model is studied with a Gaussian 
density of states of unit variance, and the spin and charge 
peaks are found to merge at U « 1.5. The fact that the 
band-width is undefined complicates comparisons with 
results in finite dimension, but one can still examine the 
ratio of U to the kinetic energy per particle, which is fi- 
nite for a Gaussian density of states. The value [//t w 4 
in two dimensions has the same ratio of U to kinetic en- 
ergy as the interaction strength at which the two peak 
structure is lost in DMFT. However, an additional diffi- 
culty in the interpretation of the DMFT results, besides 
the use of a Gaussian density of states, is the restric- 
tion of the calculations to the paramagnetic phase, and 
therefore the neglect of antiferromagnetic fiuctuations. 

It is the purpose of this paper to present a detailed 
study of the temperature dependence of the local mo- 
ment and the associated features in the specific heat for 
the half-filled two-dimensional Hubbard Hamiltonian. A 
focus of our work will be on extending the strong cou- 
pling picture of the two peak structure of C(T) to in- 
termediate and weak coupling. As we shall show, the 
connection of moment formation and moment ordering 
with the high and low temperature peaks (respectively) 
in the specific heat is modified. At the same time, we will 
describe two recently developed techniques for comput- 
ing the specific heat which hold certain advantages over 



approaches previously used. These new techniques also 
allow us to compute the entropy and free energy, quan- 
tities typically not so easy to obtain with Mopte-iJarlo. 
A fascinating conclusion of the DMFT studieal3il3 con- 
cerned the existence of a universal crossing point of the 
specific heat curves for different U. We shall show such 
a crossing occurs also in two dimensions. Finally, we will 
discuss results for dynamical quantities like the density 
of states and optical conductivity and comment on their 
consistency with the local moment and specific heat. 

II. MODEL AND METHODS 
A. The Hubbard Hamiltonian 

The two dimensional Hubbard Hamiltonian is, 

+ U Y^iui^ - i)(nu - i) - /i^("iT + "U)- (1) 

i i 

Here cl^{cj„) are creation(destruction) operators for a 
fermion of spin a on lattice site i. The kinetic energy 
term includes a sum over near neighbors (i, j) on a two- 
dimensional square lattice, and the interaction term is 
written in particle-hole symmetric form so that /i = 
corresponds to half-filling (riif +nii) = 1 for all Hamilto- 
nian parameters i, U and temperatures T. We will hence- 
forth set the hopping parameter t = 1. 

Equal time quantities of interest in this paper include 
the energy, E = (H), the specific beat C — dE/dT, the 
local moment {ml) = ((ni| — riij)^},llZl and the near neigh- 
bor spin-spin correlation function (SiSi+x). To probe 
longer range magnetic order, we evaluate the structure 
factor, 

^(Q) = ;^Ee^'^-^'"^^(("iT -"u)("JT -"ji)), (2) 

ij 

where Q = (tt, tt) is the antiferromagnetic wave vector. 

We also evaluate two dynamic quantities. The density 
of states N{lu) is given implicitly from QMC data for the 
imaginary time Green's function, 

(3) 

Likewise, the optical conductivity, axx{^), is related to 
QMC data for the imaginary time current-current corre- 
lation function, 

axAr)= {Mr)jAO)) = T^^d^^^^^' 

jx{T)^ ~itY^{4+i,a(^ia " 4aCi+xJ- (4) 
icr 
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Both N{u}) and a^xi^^) are computed using the Maxi- 
mum Entropy (ME) techniques to invert the integral 
relations. 



B. Determinant Quantum Monte Carlo 

We use determinant QMC0 to evaluate the expecta- 
tion values above. This approach treats the electron- 
electron correlations exactly, and at half-filling, where 
we focus this work, is able to produce results with very 
small statistical fluctuations at temperatures low enough 
that the ground state has been reached. The technique 
is limited to finite size lattices, and we will show ap- 
propriate scaling analyses to argue that we extract the 
thermodynamic limit. 



C. Calculation of Specific Heat 

We will evaluate the specific heat in three ways@. All 
begin by using QMC to obtain En = E{Tn) and the as- 
sociated error bars 5En at a sufficiently fine grid of Nt 
discrete temperatures r„. The first approach is straight- 
forward numerical differentiation of the energy E{Tn). 
The second utilizes a fit to the numerical data for the en- 
ergy E{Tn), and the third is an approach using the ME 
method to invert the data E{T„) to obtain a spectrum of 
excitations of the system. Th^ei^st two approaches were 
introduced relatively recently.EllEa Therefore we shall de- 
scribe them in some detail. 

In our fitting method, whose results we denote by 
Ee{TlriWe match the QMC data E^ to the functional 
form] 



M 



^e(T) = ^(0) + ^Qe- 



■/3;a 



(5) 



1=1 



by adjusting the parameters A and c/ to minimize. 



1 



Nt 

E 



{E,{Tn) ~ En))^ 
[SEnf 



(6) 



The number of parameters M is chosen to be about one- 
fourth of the number of data points to be fit. Smaller 
numbers do not allow a good fit, while larger ones overfit 
the data. We find that a range of intermediate M exists 
which gives stable and consistent results. 

Calculation of C(T) tw fitting E{T) to polynomials has 
also been used recentlyja but requires at least two sepa- 
rate functions to be used at high and low temperatures. 
An advantage of Eq. 5 is that it uses a single functional 
form over the entire T range, and has the correct low and 
high temperature limits, C{T) — > 0. 

The specific heat can also be evaluated by differentiat- 
ing an expression which relates the energy to the density 
of states of Fermi and Bose excitations in the system,c3 
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^mc(T) = - / duj uj[F{l3,uj)pF{io) + B{P,uj)pb{uj)] 

J — oo 

1 , 1 



F{P,oj) 



and differentiating to get the specific heat. 



(7) 



C„,e(T) = 
dhJ Lo\ 



dT 
dT 



+ (8) 



The integral equation for E^e (T) is inverted by using the 
ME method to obtain pf{i^) and /9b(w) from the QMC 
data for E^. We denote by E^g{T) the energy obtained 
from the resulting pf{^) and pb{^)- 

This ME approach differs in philosophy from the fitting 
approach which begins with a physically reasonable func- 
tional form Ee(T) and then minimizes the deviation 
from the numerical data. Instead, ME computes the most 
probable spectrum p{uj) given the energy data E(Tn) and 
kernals F(/3, to) and B{f3, ui), without presupposing a par- 
ticular functional form. Despite this difference, we will 
show that the results of the two techniques are very sim- 
ilar, and agree quite well with numerical differentiation. 



D. The Entropy and Free Energy 

Both the ME and fitting techniques allow the specific 
heat, entropy, and free energy to be computed by the 
standard formulae, 



CiT) 



dE{T) 



dT 

siT) = £^dr, 

F{T) = E{T) - TS{T). 



(9) 



Here E{T) = EeiT) or E„,,{T). 

In the case of the fitting technique, we can evaluate 
the sum rule. 
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TV 



TT / dT 



C{T) 
T 



M 



1=1 



lA 



21n2- 5*0, 



(10) 



which ties the high temperature entropy to the logarithm 
of the dimension of the Hilbert space. The T = entropjj 
5*0 must of course vanish in the thermodynamic limit .ci 
For the 2-d Hubbard model at C/ 7^ we find the term 
5*0 vanishes even on finite lattices, but it may be present 
in other Hamiltonians. In the present work, the sum rule 
of Eq. 10 is satisfied to a few percent. For the maximum 
entropy method, a similar check is possible by integrating 

PF- 

We now turn to the results of our simulations. 
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III. EQUAL TIME CORRELATIONS- LOCAL 
MOMENT, SPECIFIC HEAT, AND MAGNETIC 
ORDER 

A. The Local Moment 

In the introduction we reviewed the standard argu- 
ment for the expected behavior of the Hubbard model 
local moment. Early determipant QMC work for the 
2-dimensional Hubbard modelH confirmed this, as did 
subsequent investigations .c3 In Fig. 1 we show that an 
examination of with a fine temperature mesh and at 
low temperatures reveals that after reaching a plateau at 
intermediate temperatures, ml djianges value again at a 
second, low temperature, scale. E3 We will come back to 
this point in more detail later, but it is worth commenting 
immediately that while the low temperature structure in 
ml is small compared to the size of the growth at high 
temperature, it occurs over a much smaller temperature 
range, and hence contributes a large peak in the specific 
heat. 

In order to determine whether this is a finite size effect, 
in Fig. 2 we show data on a range of lattice sizes from 
4x4 to 10x10. The evidence for the existence of the low 
energy scale is robust as the lattice size is increased. We 
can also make an extrapolation to the thermodynamic 
limit assuming a correction which goes as the inverse of 
the linear system size, as spin-wave theory indicates is 
appropriate for the full structure factorE3 




T 

FIG. 2. Data for the local moment for different lattice 
sizes. The low temperature feature in remains present as 
the lattice size is increased. 



Additional insight is obtained by looking at the behav- 
ior of 771? at different values of IL as in Fig. ^. The data 
of Fig. y are replotted in Fig. Q to emphasize the uni- 
versal nature of the high temperature behavior and the 
fact that the initial increase in the local moment as tem- 
perature is decreased does indeed occur at a temperature 
scale set by U. 
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FIG. 1. The temperature dependence of the local moment 
is shown at fixed U — 4 for a 6x6 lattice. In addition to rising 
at T ~ (7, as T decreases, exhibits a second structure at 
lower temperature. 
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FIG. 3. Temperature dependence of the local moment 
for different values of the on-site repulsion U on 6x6 lattices. 
The dashed line corresponds to the zero-hopping limit, for 
U = 12. The position of the low temperature feature first 
increases in temperature as U increases, but then gradually 
falls, as emphasized later in Fig. nol 
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Scaled version of the data of Fig. 

2 _ 



The plot 

shows that for high T the local moment = 5 + ^1 that 
is, the deviation of the local moment from its noninteracting 
value ^ exhibits a universal behavior with a temperature scale 
U. 
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aligned. Localization is thereby increased with increas- 
ing temperature, giving rise to the maximum at T 0. 
This maximum at large V has also been observed in one- 
dimension. eI 

At weak coupling, we see the opposite effect. The local 
moment has instead an additional increase at low tem- 
perature. This has a natural explanation in terms of the 
formation of local magnetic order. If there is an energy 
gain with ordering, there will be an associated prefer- 
ence for large mopients. It is interesting to note that 
the DMFT resultalJ do not observe this additional mo- 
ment enhancement at weak coupling. Instead, the local 
moment always has a maximum as a function of temper- 
ature. This is, perhaps, a consequence of restricting the 
DMFT to the paramagnetic phase. 



B. The Energy and Specific Heat 

Figure ^ shows the QMC results for E{Tn) together 
with the exponential fit Ee{T). Calculation of the specific 
heat brings out the low temperature features in E{T). 
We begin our analysis of this data for the specific heat by 
looking at the data at relatively strong coupling {U = 10) 
as shown in Fig. ^. It is seen that the results for the 
Hubbard model are beautifully fit by combining the zero- 
hopping t — Q specific heat, which lies right on the high 
T Hutjbard model results, and the Heisenberg specific 
heat,E3 which similarly lies right on the low T Hubbard 
model results. The areas under both the low and high 
T peaks are precisely In 2, as expected for the high T 
loss of entropy associated with moment formation and 
then low T moment alignment. Clearly, this provides a 
good understanding of the strong coupling specific heat, 
as well as demonstrates the reliability of our approach to 
computing C{T). 
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FIG. 5. QMC values (circles) for the energy, and fittings 
provided by Eq. 5. Here U = A and the lattice is 6x6. The 
different lines (dotted, dashed, dot-dashed) show that the fit 
is stable over a range of values of the number of fitting pa- 
rameters. 



In the zero-hopping (< — 0) limit, {frij,) = 
l/{exp[— J7/ (2T)] -|- 1}, which is monotonic with T, drop- 
ping from 1 at T = to 1/2 at r ^ 00. Why is the 
maximum in {mj,) shifted from T = when the hopping 
is present? The ground state is antifcrromagnctic, with 
the effective exchange J arising from virtual hopping of 
the electrons. This virtual transfer reduces the degree 
of localization as can be seen from the values of (m^) 
at T = in Fig. |3[ In the low-lying excited states the 
deviations from the antifcrromagnctic state reduce the 
virtual hoppings since the Pauli principle forbids hop- 
ping when adjacent electron spins are ferromagnetically 
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FIG. 6. Results for the specific heat at [/ = 10 and lat- 
tice size 6x6. The circles are numerical differentiation and 
the solid line is from the exponential fitting approach. The 
dotted line is the result for C{T) at t = (that is, a single 
site). The dashed line is C{T) for the Heisenberg model with 
J = At^/U = 0.4. 
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Further confirmation of the accuracy of our C (T) cal- 
culations is evident by comparing results for the numeri- 
cal differentiation, the exponential fit, and the ME tech- 
niques as shown in Figs. ^ and ^. It is seen that the 
agreement between all three approaches is good. Figs. 
and ^ also emphasize that both the exponential fit and 
ME techniques are well suited to capturing the two en- 
ergy scales in the problem. It is also interesting to com- 
ment on the U dependence of the areas under the specific 
heat curves of Fig. ||, 0, and ||. As U decreases into the 
weak coupling regime, the low T peak has less and less 
entropy. For U = 2 the area is only about In 2/2. This 
will be discussed at greater length shortly. 
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FIG. 7. Results for the specific heat ai U — 4 and lattice 
size 6x6. The exponential fitting (full line) and ME results 
(dashed line) are in good agreement with direct numerical 
differentiation of the QMC data. Both smooth out the noise 
associated with direct numerical differentiation of the QMC 
data, though ME appears to broaden the results perhaps a 
bit too much. 



FIG. 8. Results for the specific heat at 17 = 2 and lattice 
size 6x6. The exponential fitting (full line) and ME results 
(dotted line) are both in good agreement with direct numer- 
ical differentiation of the QMC data. As at i7 = 4, ME pro- 
duces somewhat broader peaks. The dashed line is the specific 
heat for the non-interacting limit {U = 0). 



When the specific heat curves for different U are plot- 
ted together, as in Fig. ||, one sees that there is a 
nearly universal crossing at high temperature. There 
has been considerable recent discussion of this phe- 
nomenon, both its occurrence in experimental systems 
like '^He and heavy fermipa jSystems, and models like 
the Hubbard Hamiltonian.otj In the case of the Hub- 
bard model, the crossing has been argued to follow from 
the fact that the high temperature entropy is indepen- 
dent of U, ln4 = C(T, U)dT/T, which implies that 
= dC/dUdT/T. Hence dC/dU must be positive 
for some temperature ranges atid negative for others, a 
condition for crossing to occur .ta The narrowness of the 
crossing region is traced ultimately to the linear temper- 
ature dependence of the double-pccupancy, the conju- 
gate variable associated with C/Jl3 In DMFT, two cross- 
ings were observed for the Hubbard model, with the high 
temperature one being nearly universal, while the low 
temperature intersections were considerably more spread 

out. . . □ . . 

Previous studies in two dimensional exhibit crossings 
of the specific heat at = 1.6, with a crossing region 
AT* = 0.2. In Fig. ^ we confirm this result that a specific 
heat crossing occurs in two dimensions. While the cross- 
ings in the earlier studyd shift systematically with U, we 
instead see a random fiuctuation of the crossing point. 
This suggests that the width of the crossing we report 
here is dominated by statistical fiuctuations as opposed 
to possible systematic effects. We have also verified that 
the double occupancy has a linear temperature depen- 
dence at low T, especially at weak coupling, which is the 
criterion established for a universal crossing point. 
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FIG. 9. The specific lieat curves for different U values 
show a nearly universal high temperature crossing for the two 
dimensional Hubbard model, as has previously been observed 
in DMFT and in several experimental systems. 



The position of the two peaks as a function of U is 
shown in Fig. The strong coupling analysis gives 
us firm predictions for the peak positions at large U: 
First, the t = result tells us that the high T peak is 
at Thigh ~ J7/4.8. The deviation from this limit seen in 
Fig. 1C| can be ascribed to quantum fluctuations. Mean- 
while, the Heisenberg result for the low T peak is at 
Tiow = 2J/3 = 8t^/3UJB We similarly understand the 
position of Thigh at weak coupling from the U ~ anal- 
ysis: Thigh ~ i = 1 in units where t = 1. The value 
of Tiow for weak coupling is somewhat more problem- 
atic. In three dimensions, the Neel temperature which 
describes the onset of long-range niagpetic order, r,has 
a non-monotonic behavior with t/,ua first risingQ at 
small U as T/v oc exp[—2TTt/U] and subsequently faUing 
back down as T/v oc t'^/U at large U. In lower dimen- 
sions, like the 2-d case studied in this paper, T^r = 0. 
Nevertheless, in weak coupling, both the random phase 
approximatiorO and Hartree-Fock calculations give a fi- 
nite Tjv oc exp[— 271 ^yt/U] in 2~d. It is tempting in this 
case to interpret this energy scale as that of the short 
range spin fluctuations which give rise to the low tem- 
perature peak in C at weak coupling, and similarly t^/[/ 
as the corresponding energy scale at strong coupling. In- 
deed, both the increase for small U, consistent with the 
exponential form, and the subsequent decrease can be 
seen in Tiow in Fig- It might also be noted that the 
entropy in the low-T Hartree-Fock C(T) peak goes to 
as [/ ^ 0, which is also consistent with the decreasing 
entropy under the low-T QMC C (T) peak as U becomes 
small.E3 

Finite size effects are illustrated in Fig. |ll|, which shows 
the data for the specific heat, obtained by finite differ- 
ention of the energy data, on 6x6 and 10x10 lattices at 
U = 2. As is seen, the error bars as inferred from the 
scatter in the data are of the same size as any possible 
systematic effect. In determinant QMC, finite size ef- 
fects are largest at weak coupling, so this data represents 
a rather stringent test of possible lattice size dependence 
of our results for the thermodynamics. 

We now turn to the issue of the separate contributions 
of the kinetic K and potential P energies to the specific 
heat. As discussed in the introduction, we might asso- 
ciate the charge peak with the potential energy, since the 
energy U is what enforces double occupancy and reduces 
charge fluctuations. Since the energy scale J = At^ /U 
arises from virtual hopping processes, it is more natu- 
rally associated with the kinetic energy. At strong cou- 
pling this division describes the mapping of the Hubbard 
model onto the Heisenberg model, and then the speciflc 
heat of the Heisenberg model itself, and works extremely 



well quantitatively, as seen in Fig. ^. 
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FIG. 10. Position of the high T and low T peaks of the spe- 
cific heat. The dotted line is the Heisenberg limt T ~ 2J/3, 
the full line corresponds to the t = limit, T ~ (7/4.8 and 
the dash-dotted line corresponds to an RPA-like form for the 
temperature scale of the antiferromagnetic spin fluctuations, 
T - eiip[-2T\ JtJU]. 
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FIG. 11. Comparison of results for the speciflc heat at 
U = 2 and lattice sizes 6x6 and 10x10. 



At intermediate coupling it is not so natural to consider 
separately the derivatives of the kinetic and potential en- 
ergies, as these quantities mix dramatically. Indeed, the 
behavior of the local moment shown in Fig. 1 indicates 
that a,t U — 4 the potential energy in fact contributes 
both to the low and high temperature structure of C(T). 

Figure |l| shows dP/dT and dK/dT for U = 2,4, 10. 
At strong coupling U = 10, dP/dT has a high temper- 
ature maximum, while dK/dT has a low temperature 
peak, as expected. In the combined specific heat, then, 
P and K are responsible for the charge and spin peaks 



7 



respectively. Interestingly, however, even at large 
dP/dT has a significant negative dip at low T, reflect- 
ing the potential energy cost of delocalization. As we 
have remarked, this effect has previcuisly been noted in 
the 1-d Hubbard model and DMrTEM 

The positions of the contributions of dP/dT and 
dK/dT to C(T) are exchanged as U is decreased. Fi- 
nally, at L/ = 2, it is the potential energy which is re- 
sponsible for the low temperature 'spin' peak, and the 
kinetic energy for the high temperature 'charge' peak. 



0.50 




0.01 0.1 



FIG. 12. The separate temperature derivatives of tlie ki- 
netic and potential energy. At strong coupling, the tempera- 
ture dependence of the potential energy and kinetic energies 
give rise to the 'charge' and spin peaks in C respectively. This 
attribution is interchanged at weak coupling. 



The entropy S can also be obtained as the area under 
(1/T)C(T), as well as the separate kinetic and poten- 
tial contributions. Fig. ^ shows the results dX U — 2, 
[7 = 4 and U = 10. In all cases the value of S at high 
temperature equals the expected 2 In 2 to within a few 
percent. At strong coupling, U = 10, the offsetting ki- 
netic and potential contributions seen in the low-T peak 
region in Fig. ^ arc reflected as well in Fig. ^ Never- 
theless, the total entropy (solid curve) shows a shoulder 
at In 2 and then the final high-temperature value of 2 In 2 
reflecting the two peaks in Fig. |^, which correspond first 
to the "Heisenberg" disordering of the spins and then at 
higher temperature to the destruction of local moments. 
At weak coupling, the initial increase in entropy at low T 
comes from the temperature dependence of the potential 
energy which, as we have seen, is what gives rise to the 
low T peak in C{T). The area under the low T peak in 



C{T)/T is reduced from its large U value of In 2. Cer- 
tainly one origin of this decrease is the reduction of the 
local moment ml from its large U value = 1 at small 
U, as seen in Fig. ||. The entropy associated with-local 
ordering of moments scales with the moment size.E^ 
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FIG. 13. The temperature dependence of the entropy. 



C. Magnetic Correlations 

It is natural to associate the low temperature feature in 
the local moment and therefore the low T peak in C(T) 
with the onset of local antiferromagnetic correlations be- 
tween neighboring spins. To understand how these corre- 
lations develop we have calculated the spin-spin correla- 
tion function between neighboring sites {SiSj+x) and the 
magnetic structure factor 5(Q). Figure |l^ shows these 
two quantities as a function of temperature for U = 2. 
The inset shows the derivatives of the spin-spin correla- 
tion function for neighboring sites {dSiSi+x/dT) and for 
the magnetic structure factor {dS{Cl)/dT) with respect 
to the temperature: the sharp peaks form roughly at the 
same position as the specific heat has its low T peak. 
Since Tn — 0, the peak in dS{Q)/dT should not be asso- 
ciated with long ranged correlations. Indeed, the largest 
contribution to dS{Q)/dT comes from antiferromagnetic 
correlations between neighboring spins. 
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This "pseudogap" in the density of states is one of the 
central features of the 2-d Hubbard model under recent 
discussion, since it is one of the most interesting features 
of the normal state of the high tpiiperature superconduc- 
tors in the underdoped regime.c3 As the pseudogap has 
a d-wave symmetry, like the superconducting order pa- 
rameter itself, it is believed to arise as a result of short- 
range spin fluctuations which might also play a role in 
the pairing. The pseudogap's existence has a long his- 
tory of discussion, and debate, in the numerical literature 
on the Hubbard model which we shall now review, since 
concerns about possible flnite size effects in the pseudo- 
gap may bear on similar concerns in the behavior of the 
specific heat. 



FIG. 14. Spin-spin correlation between neighboring sites 
and magnetic structure factor as a function of temperature 
for a 6x6 lattice with U — 2. The insets show the derivatives 
of these two quantities with temperature. 



IV. DENSITY OF STATES AND OPTICAL 
CONDUCTIVITY 

A. The Density of States: 

While the standard interpretation of the two peak 
structure in the specific heat in terms of the freezing out 
of charge and spin fluctuations, at high and low T re- 
spectively, is consistent with our data at large [/, Fig. |l^ 
emphasized that such a picture is not as useful at weak 
coupling. 

Greater insight into the physics behind the speciflc 
heat is obtained by looking at the dynamics. Fig- 
ure |l^ shows results for the density of states N{uj) for 
U — 1, 2, 4, 6 and decreasing temperatures. At high T the 
density of states consists of a single, very broad, bump 
with maximum at w = 0. As T is lowered for the smaller 
U values (e.g., U — 2), N{uj = 0) first increases as a 
quasiparticle peak develops at a; = 0. This peak appears 
to be very similar to that found in multiband models like 
the periodic Anderson model, where it is associated with 
a Kondo resonance. As T is lowered yet further, N{uj) 
then begins to decrease as a dip begins to form in the 
center of the quasiparticle peak. 

For larger U, on the other hand (e.g., U — 6), only 
the dip develops with decreasing temperature, and so 
N{u! = 0) always decreases as T is lowered. This be- 
havior is emphasized in Fig. |l^ which shows N{uj — 0) as 
a function of temperature for U — 2,A, 6. The tempera- 
ture at which N{uj = 0) becomes small is seen in Fig. |l^ 
to increase with U, and it appears correlated with the 
position of the low T peak in C{T). Indeed, the deriva- 
tive of N{uj = 0) is maximum at the same temperature 
where C{T) has its low T peak (indicated by arrows in 
the plot). 




FIG. 15. The temperature evolution of the density of 
states at (7 = 1,2,4,6 on a 6x6 lattice. For U = 1,2, as T 
is lowered, the single broad peak first evolves into a sharper 
quasiparticle peak before a pseudogap opens. The quasipar- 
ticle peak is washed out as U increases. 
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FIG. 16. The zero frequency density of states at 
U = 2,4,6 and 6x6 lattice. The arrows indicate the location 
of the low T peaks in the specific heat. 



The dynamical data presented in Figs. p^[-fi6| do not 
make a fully compelling case for the relation between 
the low T peak in C(T) and the pseudogap, since they 
are for a single fixed lattice size. A particular issue is 
the behavior of the density of states iV(aj) of the half- 
filled Hubbard model in the thermodynamic limit. Quan- 
tum Monte Carlo results concerning this question are still 
evolving. Early simulations had the somewhat surprising 
conclusion that a pseudogap in N{uj) was present at weak 
to intermediate coupling only at T = in the thermo- 
dynamic limit. That is, while on a fixed lattice size L a 
gap in N{uj) develops at a finite temperatuce T, it would 
go away if the lattice size were increased. til Meanwhile, 
at strong coupling, the same work found the pseudogap 
persists at finite T even as the system size increases. This 
behavior was interpreted as reflecting the fact that long 
range antifcrromagnetic correlations are present only at 
T = 0, and that such long range correlations were re- 
quired for a gap in N(uj) for small U. This interpreta- 
tion was questioned, however, since one might expect the 
pseudogap to depend only on the existence of short range 
antiferromagnetic correlations. Such local order should 
form at a temperature which is independent of lattice 
size, leading to the conclusion that the pseudogap should 
be present below that temperature even on large lattices. 

If the original suggestion that the finite temperature 
pseudogap disappears at weak coupling in the thermo- 
dynamic limit were the case, it might raise similar ques- 
tions about possible finite size effects in our results for the 
low temperature structure of the magnetic moment and 
specific heat. We believe this is not a concern for three 
reasons. First, one can consider the limit of very weak 
coupling. As we already see in Fig. H, the high T peak 
in C(r) is well fit by a noninteracting calculation, and 
specifically therefore comes from the kinetic energy. If 
the local moment (potential energy) did not evolve at low 
T, then C(T) would have a single peak structure. There- 



fore, our separate finite size scaling analysis for the mo- 
ment and the specific heat support each other. Second, 
we have presented data at different finite sizes (Fig. pT| ) 
which show no evidence for the low T peak shifting with 
increased lattice size. Finally, recent work suggests that 
the pseudogap exists in the thermodynamic limji^t weak 
coupling and is not a finite size effect there.EZrc3 



B. The Optical Conductivity 

The density of states itself does not present a com- 
plete picture of the nature of the excitation gap. A more 
refined view may be obtained by looking at the optical 
conductivity, as shown in Figs. 17-19. These results show 
that the Hubbard model has a non-zero charge gap, even 
at weak to intermediate coupling where U is less than the 
bandwidth. We examine the dynamic spin susceptibility 
(not shown) and find the spin gap to vanish since the 2-d 
Hubbard model has long range magnetic order at T = 
and hence ungapped, power law, spin wave excitations. 
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FIG. 17. The frequency dependence of the optical con- 
ductivity for U — 2 and different temperatures. Inset: The 
zero frequency value as a function of temperature, indicating 
the opening of a Mott gap. 
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open the question of why the separation of the specific 
heat energy scales is maximal in intermediate dimension. 
Possibly the half-filled two-dimensional Hubbard model 
is unique due to the unusual characteristics of its nonin- 
teracting states at the Fermi energy.El However, it should 
also be recalled that the DMFT studies have been re- 
stricted to the paramagnetic phase, which may have an 
important impact on the existence of two well defined 
peaks. Finally the existence of the Nagaoka state in two 
dimensions, but not in one or infinite dimensions, indi- 
cates that there may be no reason to expect systematic 
behavior here as a function of dimension.cil 

We have emphasized that the standard nomenclature 
which identifies the high temperature peak in the specific 
heat as due to "charge" fluctuations and the low temper- 
ature peak as due to "spin" fluctuations while useful at 
strong coupling, needs to be refined. At weak coupling, 
the high temperature peak comes from the kinetic energy 
while the low temperature peak comes from the potential 
energy. By comparing with the behavior of the density of 
states, we have argued that the structure of C(T) may be 
associated with "pseudogap" formation, that is, the on- 
set of short range antiferromagnetic correlations between 
near-neighbor spins. 

A detailed understanding of the relationship of the en- 
ergy and local moment formation is desirable in a num- 
ber of contexts in using model Hamiltonians to describe 
strongly correlated materials. In particular, while mini- 
mization of the energy determines the dominant low tem- 
perature phases, the various types of behavior of the local 
moment (for example screening by conduction electrons 
in multi-band models) can also provide important clues 
concerning the suitability of diffocent models in describ- 
ing the low temperature physics c3 
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V. CONCLUSIONS 

In this paper we have examined carefully the low tem- 
perature structure of the local moment and specific heat 
of the two dimensional Hubbard model. A striking con- 
clusion of our work is that the two peak structure in 
the specific heat is preserved down to U = t = VF/8, 
where = 8t is the bandwidth. In one dimension, the 
coalescence of the spin and charge peaks appears to oc- 
cur at much larger J7, namely the one dimensional band- 
width. Meanwhile, in infinite dimension, the peaks also 
come together at an interaction strength associated with 
an average kinetic energy per particle corresponding to 
half the two dimensional bandwidth. However, we have 
pointed out that the argument that the peaks come to- 
gether which is based on a comparison of the scales At^ /U 
and U might be flawed, as the weak coup ling spin en- 
ergy scale is instead set by teyip{—2iT^Jt/U). This leaves 
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